#include "cmatrixbase.h"

/**
 * operator+= - adds given matrix to this matrix
 * @param m
 * @return this
 */
template <typename T>
CMatrixBase<T>& CMatrixBase<T>::operator += (const CMatrixBase<T> &m) {
  for (TIndex i = 0; i < this->getNumOfRows(); i++) {
    for (TIndex j = 0; j < this->getNumOfCols(); j++) {
      items[i][j] += m.getItemsPtr()[i][j];
    }
  }
  return *this;
}
//---------------------------------------------------------------------------

/**
 * operator-= - subtract given matrix from this matrix
 * @param m
 * @return this
 */
template <typename T>
CMatrixBase<T>& CMatrixBase<T>::operator -= (const CMatrixBase<T> &m) {
  for (TIndex i = 0; i < this->getNumOfRows(); i++) {
    for (TIndex j = 0; j < this->getNumOfCols(); j++) {
      items[i][j] -= m.getItemsPtr()[i][j];
    }
  }
  return *this;
}
//---------------------------------------------------------------------------

/**
 * operator*= - multiplies this matrix by given matrix
 * @param m
 * @return this
 */
template <typename T>
CMatrixBase<T>& CMatrixBase<T>::operator *= (const CMatrixBase<T> &m) {

  // Result of multiplication
  T** tmpItems = createItems(this->getNumOfRows(), m.getNumOfCols());

  // Matrix multiplication
  for (TIndex i = 0; i < this->getNumOfRows(); i++) {
    for (TIndex j = 0; j < m.getNumOfCols(); j++) {
      T tmp(0);
      for (TIndex k = 0; k < this->getNumOfCols(); k++) {
        tmp += items[i][k] * m.getItemsPtr()[k][j];
      }
      tmpItems[i][j] = tmp;
    }
  }

  // Replaces old elements with new one and return this
  reset(this->getNumOfRows(), m.getNumOfCols(), tmpItems);
  return *this;
}
//---------------------------------------------------------------------------

/**
 * operator*= - multiplies this matrix by given number
 * @param num
 * @return this
 */
template <typename T>
CMatrixBase<T>& CMatrixBase<T>::operator *= (T num) {
  for (TIndex i = 0; i < this->getNumOfRows(); i++) multiplyRow(i, num);
  return *this;
}
//---------------------------------------------------------------------------

/**
 * operator/= - divides this matrix by given number
 * @param num
 * @return this
 */
template <typename T>
CMatrixBase<T>& CMatrixBase<T>::operator /= (T num) {
  for (TIndex i = 0; i < this->getNumOfRows(); i++) divideRow(i, num);
  return *this;
}
//---------------------------------------------------------------------------

/**
 * operator== - compares elements of given matrices
 * @param m
 * @return true if all corresponding elements are equal
 */
template <typename T>
bool CMatrixBase<T>::operator == (const CMatrixBase<T> &m) const {

  // If matrices has different sizes, they are not equal
  if (this->getNumOfRows() != m.getNumOfRows() ||
      this->getNumOfCols() != m.getNumOfCols()) {
    return false;
  }

  // Iterate through elements and compare them
  for (TIndex i = 0; i < this->getNumOfRows(); i++) {
    for (TIndex j = 0; j < this->getNumOfCols(); j++) {
      if (items[i][j] != m.getItemsPtr()[i][j]) return false;
    }
  }
  return true;
}
//---------------------------------------------------------------------------

/**
 * operator!= - compare elements of given matrices
 * @param m
 * @return true if matrices are not equal
 */
template <typename T>
bool CMatrixBase<T>::operator != (const CMatrixBase<T> &m) const {
  return !(*this == m);
}
//---------------------------------------------------------------------------

/**
 * init - create identity matrix with given size
 * @param numOfRows number of rows
 * @param numOfCols number od columns
 */
template <typename T>
void CMatrixBase<T>::init(TIndex numOfRows, TIndex numOfCols) {
  setNumOfRows(numOfRows);
  setNumOfCols(numOfCols);
  items = createItems(numOfRows, numOfCols);
  clearToIdentity();
}
//---------------------------------------------------------------------------

/**
 * init - create matrix from array of elements
 * @param numOfRows
 * @param numOfCols
 * @param items
 */
template <typename T>
void CMatrixBase<T>::init(TIndex numOfRows, TIndex numOfCols, T items[]) {

  // Save dimensions
  setNumOfRows(numOfRows);
  setNumOfCols(numOfCols);

  // Create pointers to rows
  this->items = new T*[numOfRows];

  // Save elements to matrix
  for (TIndex i = 0; i < numOfRows; i++) {
    this->items[i] = new T[numOfCols];
    for (TIndex j = 0; j < numOfCols; j++) {
      this->items[i][j] = items[(i * numOfCols) + j];
    }
  }
}
//---------------------------------------------------------------------------

/**
 * init - create matrix from vector of CVectors
 * @param type type of CVectors (row/column)
 * @param v
 */
template <typename T>
void CMatrixBase<T>::init(int type, vector<CVector<T> > v) {

  // Save dimensions
  setNumOfRows(v.size());
  setNumOfCols(v[0].getSize());

  // Create pointer to rows
  items = new T*[this->getNumOfRows()];

  // Save CVectors as rows
  for (TIndex i = 0; i < this->getNumOfRows(); i++) {
    items[i] = new T[this->getNumOfCols()];
    memcpy(items[i], v[i].getItemsPtr(), this->getNumOfCols() * sizeof(T));
  }

  // If CVectors are columns, transpose matrix
  if (type == VECTOR_COL) transpose();
}
//---------------------------------------------------------------------------

/**
 * init - creates matrix from array of CVector
 * @param type type of CVectors (row/column)
 * @param numOfVectors number of CVectors
 * @param v
 */
template <typename T>
void CMatrixBase<T>::init(int type, TIndex numOfVectors, CVector<T> v[]) {

  // Save dimensions
  setNumOfRows(numOfVectors);
  setNumOfCols(v[0].getSize());

  // Create pointer to rows
  items = new T*[this->getNumOfRows()];

  // Save CVectors as rows
  for (TIndex i = 0; i < this->getNumOfRows(); i++) {
    items[i] = new T[getNumOfCols()];
    memcpy(items[i], v[i].getItemsPtr(), getNumOfCols() * sizeof(T));
  }

  // If CVectors are columns, transpose matrix
  if (type == VECTOR_COL) transpose();
}
//---------------------------------------------------------------------------

/**
 * init - creates matrix from matrix joined by vector (CVector)
 * @param m matrix
 * @param v vector
 * @param vType type of vector (row/column)
 */
template <typename T>
void CMatrixBase<T>::init(const CMatrixBase<T> &m, const CVector<T> &v, 
int vType) {

  // Temporary variables
  TIndex mr = m.getNumOfRows();
  TIndex mc = m.getNumOfCols();
  TIndex vs = v.getSize();

  // Adds new column to matrix (join matrix "m" from right)
  if (vType == this->VECTOR_COL) {

    // Number of rows and columns of resulting matrix
    setNumOfRows(mr > vs ? mr : vs);
    setNumOfCols(mc + 1);

    // Copy vector items to resized matrix, if matrix "m" has different number
    // of rows from vector "v", remaining elements of resulting matrix will be
    // zeros
    items = new T*[getNumOfRows()];
    for (TIndex i = 0; i < getNumOfRows(); i++) {

      // Create actual row of resulting matrix
      items[i] = new T[getNumOfCols()];

      // Copy actual row of matrix "m" to resulting matrix
      if (i < mr) memcpy(items[i], m.getItemsPtr()[i], mc * sizeof(T));
      else fill_n(items[i], mc, T(0));

      // Save item from actual row of given vector "v" to last column of
      // resulting matrix (row is same)
      items[i][mc] = i < vs ? v.getItemsPtr()[i] : 0;
    }
  } else {

    // Number of rows and columns of resulting matrix
    setNumOfRows(mr + 1);
    setNumOfCols(mc > vs ? mc : vs);

    // Copy elements to resulting matrix row by row, if matrix "m" has
    // different number of columns from vector "v", remaining elements of
    // resulting matrix will be zeros
    items = new T*[getNumOfRows()];
    for (TIndex i = 0; i < getNumOfRows(); i++) {

      // Create actual row of resulting matrix
      items[i] = new T[getNumOfCols()];

      // Copy actual row from given matrix "m" and vector "v" and fill empty
      // columns with zeros
      if (i == mr) {
        memcpy(items[i], v.getItemsPtr(), vs * sizeof(T));
        if (mc > vs) fill(items[i] + vs, items[i] + mc, T(0));
      } else {
        memcpy(items[i], m.getItemsPtr()[i], mc * sizeof(T));
        if (mc < vs) fill(items[i] + mc, items[i] + vs, T(0));
      }
    }
  }
}
//---------------------------------------------------------------------------

/**
 * init - creates matrix from two matrices (m1|m2)
 * @param m1 left matrix
 * @param m2 right matrix
 */
template <typename T>
void CMatrixBase<T>::init(const CMatrixBase<T> &m1, const CMatrixBase<T> &m2) {

  // Joins together two matrices "m1" and "m2" into matrix (m1|m2), if
  // matrices differs in number of rows, empty rows will be filled by zeros

  // Temporary numbers of rows and columns
  TIndex m1r = m1.getNumOfRows();
  TIndex m1c = m1.getNumOfCols();
  TIndex m2r = m2.getNumOfRows();
  TIndex m2c = m2.getNumOfCols();
    
  // Number of rows and columns of resulting matrix
  setNumOfRows(m1r > m2r ? m1r : m2r);
  setNumOfCols(m1c + m2c);

  // Copy elements
  items = new T*[getNumOfRows()];
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    items[i] = new T[getNumOfCols()];

    // Prekopiruje prvky z aktualniho radku prvni matice do odpovidajiciho
    // radku nove matice, pokud aktualni radek presahuje pocet radku prvni
    // matice, vyplni radek nulami

    // Copy elements from actual row of first matrix "m1" if index of actual
    // row overlaps number of rows in "m1", fill row with zeros
    if (i < m1r) memcpy(items[i], m1.getItemsPtr()[i], m1c * sizeof(T));
    else fill_n(items[i], m1c, T(0));

    // Copy elements from actual row of second matrix "m2" if index of actual
    // row overlaps number of rows in "m2", fill row with zeros
    if (i < m2r) memcpy(items[i] + m1c, m2.getItemsPtr()[i], m2c * sizeof(T));
    else fill_n(items[i] + m1c, m2c, T(0));
  }
}
//---------------------------------------------------------------------------

/**
 * getRow - gets row from matrix
 * @param index index of row
 * @return row at position "index"
 */
template <typename T>
CVector<T> CMatrixBase<T>::getRow(TIndex index) const {
  if (index >= getNumOfRows()) {
    throw CIndexOutOfBoundsException(getClassName(), "getRow");
  }
  CVector<T> result(getNumOfCols());
  for (TIndex i = 0; i < getNumOfCols(); i++) result.set(i, items[index][i]);
  return result;
}
//---------------------------------------------------------------------------

/**
 * getCol - gets column from matrix
 * @param index index of column
 * @return column at position "index"
 */
template <typename T>
CVector<T> CMatrixBase<T>::getCol(TIndex index) const {
  if (index >= getNumOfCols()) {
    throw CIndexOutOfBoundsException(getClassName(), "getCol");
  }
  CVector<T> result(getNumOfRows());
  for (TIndex i = 0; i < getNumOfRows(); i++) result.set(i, items[i][index]);
  return result;
}
//---------------------------------------------------------------------------

/**
 * get - gets value of element in given row and column
 * @param row row index
 * @param col column index
 * @return element at position "[row;col]"
 */
template <typename T>
T CMatrixBase<T>::get(TIndex row, TIndex col) const {
  if (row >= getNumOfRows() || col >= getNumOfCols()) {
    throw CIndexOutOfBoundsException(getClassName(), "get");
  }
  return items[row][col];
}
//---------------------------------------------------------------------------

/**
 * set - sets value for element in given row and column
 * @param row row index
 * @param col column index
 * @param value
 */
template <typename T>
void CMatrixBase<T>::set(TIndex row, TIndex col, T value) {
  if (row >= getNumOfRows() || col >= getNumOfCols()) {
    throw CIndexOutOfBoundsException(getClassName(), "set");
  }
  items[row][col] = value;
}
//---------------------------------------------------------------------------

/**
 * swapRows - swaps two rows of matrix
 * @param index1
 * @param index2
 */
template <typename T>
void CMatrixBase<T>::swapRows(TIndex index1, TIndex index2) {
  if (index1 >= getNumOfRows() || index2 >= getNumOfRows()) {
    throw CIndexOutOfBoundsException(getClassName(), "swapRows");
  }
  T* tmp = items[index1];
  items[index1] = items[index2];
  items[index2] = tmp;
}
//---------------------------------------------------------------------------

/**
 * multiplyRow - multiply one row of matrix by given number
 * @param index index of row
 * @param num factor
 */
template <typename T>
void CMatrixBase<T>::multiplyRow(TIndex index, T num) {
  if (index >= getNumOfRows()) {
    throw CIndexOutOfBoundsException(getClassName(), "multiplyRow");
  }
  for (TIndex i = 0; i < getNumOfCols(); i++) items[index][i] *= num;
}
//---------------------------------------------------------------------------

/**
 * divideRow - divide one row of matrix by given number
 * @param index index of row
 * @param num divisor
 */
template <typename T>
void CMatrixBase<T>::divideRow(TIndex index, T num) {
  if (index >= getNumOfRows()) {
    throw CIndexOutOfBoundsException(getClassName(), "divideRow");
  }
  for (TIndex i = 0; i < getNumOfCols(); i++) items[index][i] /= num;
}
//---------------------------------------------------------------------------

/**
 * clearToIdentity - replaces matrix elements with elements of identity matrix
 */
template <typename T>
void CMatrixBase<T>::clearToIdentity() {
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    for (TIndex j = 0; j < getNumOfCols(); j++) items[i][j] = i == j ? 1 : 0;
  }
}
//---------------------------------------------------------------------------

/**
 * clearToZero - replaces matrix elements with zeros
 */
template <typename T>
void CMatrixBase<T>::clearToZero() {
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    fill_n(items[i], getNumOfCols(), T(0));
  }
}
//---------------------------------------------------------------------------

/**
 * getItemsPtr - returns pointer to matrix elements
 * @return pointer to matrix elements
 */
template <typename T>
T** CMatrixBase<T>::getItemsPtr() const {
  return items;
}
//---------------------------------------------------------------------------

/**
 * isSymmetric - checks whether matrix is symmetric
 * @return true if matrix is symmetric
 */
template <typename T>
bool CMatrixBase<T>::isSymmetric() const {

  // If matrix is not square, it is not symmetric, return
  if (!isSquare()) return false;

  // Test whether A^T == A
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    for (TIndex j = 0; j < i; j++) {
      if (items[i][j] != items[j][i]) return false;
    }
  }
  return true;
}
//---------------------------------------------------------------------------

/**
 * isSquare - checks whether matrix is square
 * @return true if matrix is square
 */
template <typename T>
bool CMatrixBase<T>::isSquare() const {
  return getNumOfRows() == getNumOfCols();
}
//---------------------------------------------------------------------------

/**
 * isUpperTriangular - checks whether matrix is upper triangular
 * @return true if matrix is upper triangular
 */
template <typename T>
bool CMatrixBase<T>::isUpperTriangular() const {
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    for (TIndex j = 0; j < i; j++) {
      if (!CMath<T>::equals(0, items[i][j])) {
        return false;
      }
    }
  }
  return true;
}
//---------------------------------------------------------------------------

/**
 * doLUDecomposition - performs LU decomposition
 * @param l resulting lower triangular matrix
 * @param u resulting upper triangular matrix
 */
template <typename T>
void CMatrixBase<T>::doLUDecomposition(CMatrixBase<T> &l,
    CMatrixBase<T> &u) const {

  // If matrix isn't square, exception
  if (!isSquare()) {
    throw CMatrixNotSquareException(getClassName(), "doLUDecomposition");
  }
  
  // Initialize resulting matrices
  u.free();
  l.free();
  u.init(getNumOfRows(), getNumOfCols());
  l.init(getNumOfRows(), getNumOfCols());
  u.clearToZero();
  
  // Decompose matrix - Crout algorithm
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    for (TIndex j = 0; j < getNumOfCols(); j++) {            
      T sum(0);      
      if (i > j) {
        for (TIndex k = 0; k < j; k++) {
          sum += u.getItemsPtr()[k][j] * l.getItemsPtr()[i][k];
        }
        l.set(i, j, (items[i][j] - sum) / u.getItemsPtr()[j][j]);
      } else {
        for (TIndex k = 0; k < i; k++) {
          sum += u.getItemsPtr()[k][j] * l.getItemsPtr()[i][k];
        }
        u.set(i, j, items[i][j] - sum);
      }
    }
  }
}
//---------------------------------------------------------------------------

/**
 * doGaussJordanElimination - does full Gauss-Jordan elimination
 * @return true if elimination was successfull
 */
template <typename T>
bool CMatrixBase<T>::doGaussJordanElimination() {
  doGaussElimination();
  return toLowerTriangular();
}
//---------------------------------------------------------------------------

/**
 * toLowerTriangular - transforms (eliminate) matrix to lower triangular form
 * @return true if transformation was successfull
 */
template <typename T>
bool CMatrixBase<T>::toLowerTriangular() {

  // Calculate number of diagonal elements, if is zero, return
  TIndex i = getNumOfRows() < getNumOfCols() ? getNumOfRows() : getNumOfCols();
  if (!i) return false;

  // Go through diagonal elements (from back)
  do {
    TIndex j = --i;

    // If actual element is zero, return
    if (!items[i][i]) return false;

    // Eliminate all elements in column above actual element
    do {
      T tmp = -items[--j][i];
      for (TIndex k = i; k < getNumOfCols(); k++) {
        items[j][k] += tmp * items[i][k];
      }
    } while (j);
  } while (i > 1);

  return true;
}
//---------------------------------------------------------------------------

/**
 * doGaussElimination - transform (eliminate) matrix to upper triangular form
 * @return number of LD (linear dependent) columns
 */
template <typename T>
TIndex CMatrixBase<T>::doGaussElimination() {

  // Iterate through main diagonal
  TIndex i = 0, j = 0;
  while (i < getNumOfRows() && j < getNumOfCols()) {
    
    // If actual element is zero after pivoting => actual column is LD, move
    // to neighbouring right element
    doPivoting(i, j);
    if (items[i][j] == 0) {
      j++;
      continue;
    }

    // Eliminate column under actual element and move to next element
    elimColUnder(i++, j++);
  }

  // Return number of LD columns (distance between j and i plus number of
  // LD columns after upper triangular part of matrix)
  return (j - i) + (getNumOfCols() - j);
}
//---------------------------------------------------------------------------

/**
 * elimColUnder - eliminate column below element in given row and column
 * @param row
 * @param col
 * @return value which was at position "row;col" before elimination
 */
template <typename T>
T CMatrixBase<T>::elimColUnder(TIndex row, TIndex col) {

  // Test indices
  if (row >= getNumOfRows() || col >= getNumOfCols()) {
    throw CIndexOutOfBoundsException(getClassName(), "elimColUnder");
  }

  // Temporary variable
  T tmpFactor = 1;

  // If given element differs from one, divide row right from given element
  // by given element
  if (items[row][col] != 1) {
    tmpFactor = items[row][col];
    for (TIndex j = col; j < getNumOfCols(); j++) items[row][j] /= tmpFactor;
  }

  // Eliminate column below given element
  for (TIndex j = row + 1; j < getNumOfRows(); j++) {
    T tmp = -items[j][col];
    for (TIndex k = col; k < getNumOfCols(); k++) {
      items[j][k] += tmp * items[row][k];
    }
  }

  // Return original value at [row;col]
  return tmpFactor;
}
//---------------------------------------------------------------------------

/**
 * elimCol - eliminate column up and below element in given row and column
 * @param row
 * @param col
 * @return value which was at position "row;col" before elimination
 */
template <typename T>
T CMatrixBase<T>::elimCol(TIndex row, TIndex col) {

  // Test indices
  if (row >= getNumOfRows() || col >= getNumOfCols()) {
    throw CIndexOutOfBoundsException(getClassName(), "elimCol");
  }

  // Temporary value
  T tmpFactor = 1;

  // If given element differs from one divide actual row by given element
  if (items[row][col] != 1) {
    tmpFactor = items[row][col];
    for (TIndex j = 0; j < getNumOfCols(); j++) items[row][j] /= tmpFactor;
  }

  // Eliminate all elements up and below given element
  for (TIndex j = 0; j < getNumOfRows(); j++) {
    if (j == row) continue;
    T tmp = -items[j][col];
    for (TIndex k = 0; k < getNumOfCols(); k++) {
      items[j][k] += tmp * items[row][k];
    }
  }

  // Return original value at [row;col]
  return tmpFactor;
}
//---------------------------------------------------------------------------

/**
 * doPivoting - removes zero from main diagonal by swapping matrix rows
 * @param row row of pivot
 * @param col column of pivot
 * @return true if rows were swapped
 */
template <typename T>
bool CMatrixBase<T>::doPivoting(TIndex row, TIndex col) {
  
  // Finds maximal element in given column and under given row
  T max(items[row][col]);
  TIndex maxRow(row);
  for (TIndex i = row; i < getNumOfRows(); i++) {
    if (CMath<T>::abs(items[i][col]) > CMath<T>::abs(max)) {
      max = items[i][col];
      maxRow = i;
    }
  }
  
  // Swap actual row with found row (with maximal element)
  if (row != maxRow) {
    swapRows(row, maxRow);
    return true;
  }
  return false;
}
//---------------------------------------------------------------------------

/**
 * transpose - transpose matrix
 */
template <typename T>
void CMatrixBase<T>::transpose() {

  // Temporary elements
  T** tmpItems = createItems(getNumOfCols(), getNumOfRows());

  // Swap rows and columns
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    for (TIndex j = 0; j < getNumOfCols(); j++) tmpItems[j][i] = items[i][j];
  }

  // Delete matrix elements and assign swapped elements to it
  reset(getNumOfCols(), getNumOfRows(), tmpItems);
}
//---------------------------------------------------------------------------

/**
 * createItems - allocates memory for matrix elements
 * @param numOfRows
 * @param numOfCols
 * @return prvky matice
 */
template <typename T>
T** CMatrixBase<T>::createItems(TIndex numOfRows, TIndex numOfCols) {
  T** tmpItems = new T*[numOfRows];
  for (TIndex i = 0; i < numOfRows; i++) tmpItems[i] = new T[numOfCols];
  return tmpItems;
}
//---------------------------------------------------------------------------

/**
 * copyFrom - replace this matrix with copy of given matrix (deep copy)
 * @param m
 */
template <typename T>
void CMatrixBase<T>::copyFrom(const CMatrixBase<T> &m) {
  
  // Save dimensions
  setNumOfRows(m.getNumOfRows());
  setNumOfCols(m.getNumOfCols());
  
  // Copy elements
  items = new T*[getNumOfRows()];
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    items[i] = new T[getNumOfCols()];
    memcpy(items[i], m.getItemsPtr()[i], getNumOfCols() * sizeof(T));
  }
}
//---------------------------------------------------------------------------

/**
 * toString - saves matrix elements to string
 * @return string
 */
template <typename T>
string CMatrixBase<T>::toString() const {

  // Output stream
  ostringstream stream;

  // Iterate through elements and print them to output stream
  for (TIndex i = 0; i < getNumOfRows(); i++) {
    for (TIndex j = 0; j < getNumOfCols(); j++) {
      stream << items[i][j] << " ";
    }
    stream << endl;
  }
  return stream.str();
}
//---------------------------------------------------------------------------

/**
 * mul - multiplies given matrix by given vector
 * @param m matrix
 * @param v vector
 * @return product (a vector)
 */
template <typename T>
CVector<T> CMatrixBase<T>::mul(const CMatrixBase<T> &m, const CVector<T> &v) {

  // Resulting product
  CVector<T> result(v.getSize());

  // Multiplication
  for (TIndex i = 0; i < m.getNumOfRows(); i++) {
    T tmp = 0;
    for (TIndex j = 0; j < v.getSize(); j++) {
      tmp += m.getItemsPtr()[i][j] * v.getItemsPtr()[j];
    }
    result.set(i, tmp);
  }
  return result;
}
//---------------------------------------------------------------------------

/**
 * mul - multiplies given vector by given matrix
 * @param v vector
 * @param m matrix
 * @return product (a vector)
 */
template <typename T>
CVector<T> CMatrixBase<T>::mul(const CVector<T> &v, const CMatrixBase<T> &m) {

  // Resulting product
  CVector<T> result(v.getSize());

  // Multiplication
  for (TIndex i = 0; i < m.getNumOfCols(); i++) {
    T tmp = 0;
    for (TIndex j = 0; j < v.getSize(); j++) {
      tmp += m.getItemsPtr()[j][i] * v.getItemsPtr()[j];
    }
    result.set(i, tmp);
  }
  return result;
}
//---------------------------------------------------------------------------

/**
 * reset - resets matrix
 * @param numOfRows
 * @param numOfCols
 * @param items
 */
template <typename T>
void CMatrixBase<T>::reset(TIndex numOfRows, TIndex numOfCols, T** items) {
  free();
  setNumOfRows(numOfRows);
  setNumOfCols(numOfCols);
  this->items = items;
}
//---------------------------------------------------------------------------

/**
 * free - deletes matrix from memory
 */
template <typename T>
void CMatrixBase<T>::free() {
  if (items != NULL) {
    for (TIndex i = 0; i < getNumOfRows(); i++) {
      if (items[i] != NULL) delete [] items[i];
    }
    delete [] items;
    items = NULL;
  }
  setNumOfRows(0);
  setNumOfCols(0);
}
//---------------------------------------------------------------------------

template class CMatrixBase<char>;
template class CMatrixBase<short>;
template class CMatrixBase<int>;
template class CMatrixBase<long>;
template class CMatrixBase<float>;
template class CMatrixBase<double>;
template class CMatrixBase<CBigInt>;
template class CMatrixBase<CFraction>;
